﻿using Bling.DSL;
using Bling.Util;
using System.Collections.Generic;
using System;
using Bling.Core;

namespace Bling.NewOps {
  namespace Conversions {
    using Bling.NewVecs;
    public static class Extensions {
      public static Expr<IVec<bool, D>> GreaterThan<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return CompareOps<T, D>.GreaterThan.Make(E0, E1);
      }
      public static Expr<IVec<bool, D>> LessThan<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return CompareOps<T, D>.LessThan.Make(E0, E1);
      }
      public static Expr<IVec<bool, D>> GreaterThanEqual<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return CompareOps<T, D>.GreaterThanEqual.Make(E0, E1);
      }
      public static Expr<IVec<bool, D>> LessThanEqual<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return CompareOps<T, D>.LessThanEqual.Make(E0, E1);
      }
      public static Expr<IVec<T, D>> Condition<T, D>(this Expr<IVec<bool, D>> Test, Expr<IVec<T, D>> IfTrue, Expr<IVec<T, D>> IfFalse) where D : Dim<D> {
        throw new NotImplementedException();
      }
      public static Expr<IVec<bool, D>> Equals0<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return EqualityOps<T, D>.Equals.Make(E0, E1);
      }
      public static Expr<IVec<bool, D>> NotEquals0<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return EqualityOps<T, D>.NotEquals.Make(E0, E1);
      }
      public static Expr<IVec<bool, D>> And<D>(this Expr<IVec<bool, D>> E0, Expr<IVec<bool, D>> E1) where D : Dim<D> {
        return BoolOps<D>.And.Make(E0, E1);
      }
      public static Expr<IVec<bool, D>> Or<D>(this Expr<IVec<bool, D>> E0, Expr<IVec<bool, D>> E1) where D : Dim<D> {
        return BoolOps<D>.Or.Make(E0, E1);
      }
      public static Expr<IVec<bool, D>> ExclusiveOr<D>(this Expr<IVec<bool, D>> E0, Expr<IVec<bool, D>> E1) where D : Dim<D> {
        return BoolOps<D>.ExclusiveOr.Make(E0, E1);
      }
      public static Expr<IVec<bool, D>> Not<D>(this Expr<IVec<bool, D>> E0) where D : Dim<D> {
        return BoolOps<D>.Not.Make(E0);
      }
      public static Expr<IVec<bool, D1>> Any<D>(this Expr<IVec<bool, D>> E0) where D : Dim<D> {
        return BoolOps<D>.Any.Make(E0);
      }
      public static Expr<IVec<bool, D1>> All<D>(this Expr<IVec<bool, D>> E0) where D : Dim<D> {
        return BoolOps<D>.All.Make(E0);
      }
      public static Expr<IVec<T, D>> Bl<T, D>(this T t) where D : Dim<D> {
        var c = new Constant<T>(t);
        var Es = new Expr<T>[Dim<D>.Count];
        for (int i = 0; i < Es.Length; i++) Es[i] = c;
        return Composite<T, D>.Make(Es);
      }
      public static Expr<IVec<T, D1>> Repeat<T>(this T t) {
        return Bl<T, D1>(t);
      }
      public static Expr<IVec<T, D>> Bl<T, D>(this Expr<T> t) where D : Dim<D> {
        if (t is ProxyExpr<T>) // push the proxy down!
          return ((ProxyExpr<T>)t).Make0<IVec<T, D>>(e => e.Bl<T, D>());
        var Es = new Expr<T>[Dim<D>.Count];
        for (int i = 0; i < Es.Length; i++) Es[i] = t;
        return Composite<T, D>.Make(Es);
      }
      public static Expr<IVec<T, D1>> Bl0<T>(this Expr<T> t) {
        return Bl<T, D1>(t);
      }
      public static Expr<IVec<T, D>> Bl<T,D>(this Expr<IVec<T,D1>> t) where D : Dim<D> {
        return t.Access(0).Bl<T, D>();
      }
      public static Expr<T> Access<T, D>(this Expr<IVec<T, D>> E0, Expr<int> n) where D : Dim<D> {
        if (n is Constant<int>) return Access(E0, ((Constant<int>)n).Value);
        else if (Dim<D>.Count == 1) return Access(E0, 0);
        if (E0 is ProxyExpr<IVec<T, D>>) {
          return ((ProxyExpr<IVec<T, D>>)E0).Make0(d => d.Access(n));
        }
        return new AccessExpr<T, D>(E0, n);
      }
      public static Expr<T> Access<T, D>(this Expr<IVec<T, D>> E0, int n) where D : Dim<D> {
        if (E0 is ProxyExpr<IVec<T, D>>) {
          return ((ProxyExpr<IVec<T, D>>)E0).Make0(d => d.Access(n));
        }
        (n < Dim<D>.Count).Assert();
        if (E0 is Composite<T, D>) {
          return (((Composite<T, D>)E0)[n]);
        } else return new AccessExpr<T, D>(E0, n);
      }
      public static Expr<IVec<T, D>> Add<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return NumericOps<T, D>.Add.Make(E0, E1);
      }
      public static Expr<IVec<T, D>> Subtract<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return NumericOps<T, D>.Subtract.Make(E0, E1);
      }
      public static Expr<IVec<T, D>> Negate<T, D>(this Expr<IVec<T, D>> E0) where D : Dim<D> {
        return NumericOps<T, D>.Negate.Make(E0);
      }
      public static Expr<IVec<T, D>> Abs<T, D>(this Expr<IVec<T, D>> E0) where D : Dim<D> {
        return NumericOps<T, D>.Abs.Make(E0);
      }
      public static Expr<IVec<T, D>> Multiply<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return NumericOps<T, D>.Multiply.Make(E0, E1);
      }
      public static Expr<IVec<T, D>> Divide<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return NumericOps<T, D>.Divide.Make(E0, E1);
      }
      public static Expr<IVec<T, D>> Recip<T, D>(this Expr<IVec<T, D>> E0) where D : Dim<D> {
        return NumericOps<T, D>.Recip.Make(E0);
      }
      public static Expr<IVec<T, D>> Mod<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return NumericOps<T, D>.Mod.Make(E0, E1);
      }
      public static Expr<IVec<T, D>> Min<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return NumericOps<T, D>.Min.Make(E0, E1);
      }
      public static Expr<IVec<T, D>> Max<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1) where D : Dim<D> {
        return NumericOps<T, D>.Max.Make(E0, E1);
      }
      public static Expr<IVec<T, D>> Clamp<T, D>(this Expr<IVec<T, D>> E0, Expr<IVec<T, D>> E1, Expr<IVec<T, D>> E2) where D : Dim<D> {
        return NumericOps<T, D>.Clamp.Make(E0, E1, E2);
      }
      public static Expr<IVec<double, D>> Saturate<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Saturate.Make(E0);
      }
      public static Expr<IVec<double, D>> Floor<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Floor.Make(E0);
      }
      public static Expr<IVec<double, D>> Truncate<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Truncate.Make(E0);
      }
      public static Expr<IVec<double, D>> Fmod<D>(this Expr<IVec<double, D>> E0, Expr<IVec<double, D>> E1) where D : Dim<D> {
        return DoubleOps<D>.Fmod.Make(E0, E1);
      }
      public static Expr<IVec<double, D>> FrExp<D>(this Expr<IVec<double, D>> E0, Expr<IVec<double, D>> E1) where D : Dim<D> {
        return DoubleOps<D>.FrExp.Make(E0, E1);
      }
      public static Expr<IVec<double, D>> Ceiling<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Ceiling.Make(E0);
      }
      public static Expr<IVec<double, D>> Round<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Round.Make(E0);
      }
      public static Expr<IVec<double, D>> Frac<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Frac.Make(E0);
      }
      public static Expr<IVec<double, D>> Square<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Square.Make(E0);
      }
      public static Expr<IVec<double, D>> Sqrt<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Sqrt.Make(E0);
      }
      public static Expr<IVec<double, D>> Rsqrt<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Rsqrt.Make(E0);
      }
      public static Expr<IVec<double, D1>> Dot<D>(this Expr<IVec<double, D>> E0, Expr<IVec<double, D>> E1) where D : Dim<D> {
        return DoubleOps<D>.Dot.Make(E0, E1);
      }
      public static Expr<IVec<double, D1>> LengthSquared<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.LengthSquared.Make(E0);
      }
      public static Expr<IVec<double, D1>> Length<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Length.Make(E0);
      }
      public static Expr<IVec<double, D>> Normalize<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Normalize.Make(E0);
      }
      public static Expr<IVec<double, D>> Sign<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Sign.Make(E0);
      }
      public static Expr<IVec<double, D>> Faceforward<D>(this Expr<IVec<double, D>> E0, Expr<IVec<double, D>> E1, Expr<IVec<double, D>> E2) where D : Dim<D> {
        return DoubleOps<D>.Faceforward.Make(E0, E1, E2);
      }
      public static Expr<IVec<double, D>> Reflect<D>(this Expr<IVec<double, D>> E0, Expr<IVec<double, D>> E1) where D : Dim<D> {
        return DoubleOps<D>.Reflect.Make(E0, E1);
      }
      public static Expr<IVec<double, D>> Refract<D>(this Expr<IVec<double, D>> E0, Expr<IVec<double, D>> E1, Expr<IVec<double, D1>> E2) where D : Dim<D> {
        return DoubleOps<D>.Refract.Make(E0, E1, E2);
      }

      public static Expr<IVec<double, D>> Ddx<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Ddx.Make(E0);
      }
      public static Expr<IVec<double, D>> Ddy<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Ddy.Make(E0);
      }
      public static Expr<IVec<double, D>> FWidth<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.FWidth.Make(E0);
      }
      public static Expr<IVec<double, D>> Exp<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Exp.Make(E0);
      }
      public static Expr<IVec<double, D>> Exp2<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Exp2.Make(E0);
      }
      public static Expr<IVec<double, D>> Exp10<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Exp10.Make(E0);
      }
      public static Expr<IVec<double, D>> Pow<D>(this Expr<IVec<double, D>> E0, Expr<IVec<double, D>> E1) where D : Dim<D> {
        return DoubleOps<D>.Pow.Make(E0, E1);
      }
      public static Expr<IVec<double, D>> Log<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Log.Make(E0);
      }
      public static Expr<IVec<double, D>> Log2<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Log2.Make(E0);
      }
      public static Expr<IVec<double, D>> Log10<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.Log10.Make(E0);
      }
      public static Expr<IVec<double, D>> LogN<D>(this Expr<IVec<double, D>> E0, Expr<IVec<double, D>> E1) where D : Dim<D> {
        return DoubleOps<D>.LogN.Make(E0, E1);
      }

      public static Expr<IVec<bool, D1>> IsNaN<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.IsNaN.Make(E0);
      }
      public static Expr<IVec<bool, D1>> IsInfinity<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.IsInfinity.Make(E0);
      }
      public static Expr<IVec<bool, D1>> IsFinite<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.IsFinite.Make(E0);
      }
      public static Expr<IVec<bool, D1>> IsPositiveInfinity<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.IsPositiveInfinity.Make(E0);
      }
      public static Expr<IVec<bool, D1>> IsNegativeInfinity<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return DoubleOps<D>.IsNegativeInfinity.Make(E0);
      }
      public static Expr<IVec<double, D>> Sin<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Sin.Make(E0);
      }
      public static Expr<IVec<double, D>> Asin<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Asin.Make(E0);
      }
      public static Expr<IVec<double, D>> Sinh<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Sinh.Make(E0);
      }
      public static Expr<IVec<double, D>> Cos<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Cos.Make(E0);
      }
      public static Expr<IVec<double, D>> Acos<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Acos.Make(E0);
      }
      public static Expr<IVec<double, D>> Cosh<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Cosh.Make(E0);
      }
      public static Expr<IVec<double, D>> Tan<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Tan.Make(E0);
      }
      public static Expr<IVec<double, D>> Atan<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Atan.Make(E0);
      }
      public static Expr<IVec<double, D>> Atan2<D>(this Expr<IVec<double, D>> E0, Expr<IVec<double, D>> E1) where D : Dim<D> {
        return TrigOps<D>.Atan2.Make(E0, E1);
      } 
      public static Expr<IVec<double, D>> Radians<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Radians.Make(E0);
      }
      public static Expr<IVec<double, D>> Degrees<D>(this Expr<IVec<double, D>> E0) where D : Dim<D> {
        return TrigOps<D>.Degrees.Make(E0);
      }

      public static Expr<IVec<T, D>> Transpose<T,D>(this Expr<IVec<T, D>> E0)
        where D : Dim2D<D> {
        return MatrixOps<T, D>.Transpose.Make(E0);
      }
      public static Expr<T> Access<T, D>(this Expr<IVec<T, D>> E0, int m, int n)
        where D : Dim2D<D> {
        (m < Dim2D<D>.MCount).Assert();
        (n < Dim2D<D>.NCount).Assert();
        if (E0 is ProxyExpr<IVec<T, D>>) return ((ProxyExpr<IVec<T, D>>)E0).Make0(E1 => E1.Access<T, D>(m, n));
        return new Access2DExpr<T, D>(E0, m, n);
      }
      public static Expr<T> Access<T, D>(this Expr<IVec<T, D>> E0, Expr<int> m, Expr<int> n)
        where D : Dim2D<D> {
        if (E0 is ProxyExpr<IVec<T, D>>) return ((ProxyExpr<IVec<T, D>>)E0).Make0(E1 => E1.Access<T, D>(m, n));
        return new Access2DExpr<T, D>(E0, m, n);
      }
      public static T MaybeNull<T>(this bool b, T t) { return b ? t : default(T); }
    }
  }

}
namespace Bling.NewOps {
  using Bling.NewVecs;
  using Bling.NewOps.Conversions;
  public partial class BasicOps<T> {
  }
  public partial class EqualityOps<T, D> where D : Dim<D> {
    public class BinaryOp : Binary2Op<T, T, bool, D>, IPrimitiveDefined {
      public BinaryOp() {
        Defined<T, D>.Add(this);
      }
    }
    public static new readonly BinaryOp Equals = new BinaryOp() {
      Id = "==",
      MonoOp = EqualityOps<T, D1>.Equals,
      Op = EqualityOps<object,D1>.Equals,
    };
    public static readonly BinaryOp NotEquals = new BinaryOp() {
      Id = "!=",
      CustomExpand = (arg0, arg1) => arg0.Equals0(arg1).Not(),
      MonoOp = EqualityOps<T, D1>.NotEquals,
      Op = EqualityOps<object, D1>.NotEquals,
    };
  }
  public partial class CompareOps<T, D> where D : Dim<D> {
    public class BinaryOp : Binary2Op<T, T, bool, D>, IPrimitiveDefined {
      public BinaryOp() {
        Defined<T, D>.Add(this);
      }
    }
    public static readonly BinaryOp GreaterThan = new BinaryOp() {
      Id = ">",
      MonoOp = CompareOps<T, D1>.GreaterThan,
      Op = CompareOps<object, D1>.GreaterThan,
    };
    public static readonly BinaryOp LessThan = new BinaryOp() {
      Id = "<",
      MonoOp = CompareOps<T, D1>.LessThan,
      Op = CompareOps<object, D1>.LessThan,
    };
    public static readonly BinaryOp GreaterThanEqual = new BinaryOp() {
      Id = ">=",
      MonoOp = CompareOps<T, D1>.GreaterThanEqual,
      Op = CompareOps<object, D1>.GreaterThanEqual,
    };
    public static readonly BinaryOp LessThanEqual = new BinaryOp() {
      Id = "<=",
      MonoOp = CompareOps<T, D1>.LessThanEqual,
      Op = CompareOps<object, D1>.LessThanEqual,
    };
  }
  public partial class BoolOps<D> where D : Dim<D> {
    public class BinaryOp : Binary2Op<bool, bool, bool, D>, IPrimitiveDefined {
      public BinaryOp() {
        Defined<bool, D>.Add(this);
      }
    }
    public class UnaryOp : UnaryOp<bool, bool, D>, IPrimitiveDefined {
      public UnaryOp() {
        Defined<bool, D>.Add(this);
      }
    }
    public class Unary2Op : UnaryOp<IVec<bool, D>, IVec<bool, D1>>, IShaderDefined {
      public Unary2Op() {
        Defined<bool, D>.Add(this);
      }
    }
    public static readonly BinaryOp And = new BinaryOp() {
      Id = "&",
      MonoOp = BoolOps<D1>.And,
      Op = BoolOps<D1>.And,
    };
    public static readonly BinaryOp Or = new BinaryOp() {
      Id = "|",
      MonoOp = BoolOps<D1>.Or,
      Op = BoolOps<D1>.Or,
    };
    public static readonly BinaryOp ExclusiveOr = new BinaryOp() {
      Id = "^",
      MonoOp = BoolOps<D1>.ExclusiveOr,
      Op = BoolOps<D1>.ExclusiveOr,
    };
    public static readonly UnaryOp Not = new UnaryOp() {
      Id = "!",
      MonoOp = BoolOps<D1>.Not,
      Op = BoolOps<D1>.Not,
    };
    public static readonly Unary2Op Any = new Unary2Op() {
      Id = "Any",
      CustomExpand = (arg0) => {
        var result = false.Repeat<bool>();
        for (int i = 0; i < Dim<D>.Count; i++) {
          result = result.Or(arg0.Access(i).Bl0());
        }
        return result;
      },
      Op = BoolOps<D1>.Any,
    };
    public static readonly Unary2Op All = new Unary2Op() {
      Id = "All",
      CustomExpand = (arg0) => {
        var result = true.Repeat<bool>();
        for (int i = 0; i < Dim<D>.Count; i++) {
          result = result.And(arg0.Access(i).Bl0());
        }
        return result;
      },
      Op = BoolOps<D1>.All,
    };
  }
  public static partial class NumericOps<T, D> where D : Dim<D> {
    public class UnaryOp : UnaryOp<T, T, D> { }
    public class BinaryOp : Binary2Op<T, T, T, D>, IPrimitiveDefined {
      public BinaryOp() {
        Defined<bool, D>.Add(this);
      }
    }
    public class TrinaryOp : TrinaryOp<T, T, T, T, D>, IShaderDefined {
      public TrinaryOp() {
        Defined<bool, D>.Add(this);
      }
    }
    public class UnaryPOp : UnaryOp, IPrimitiveDefined, IMathDefined {
      public UnaryPOp() {
        Defined<bool, D>.Add(this);
      }
    
    }
    public class UnaryQOp : UnaryOp, IMathDefined, IShaderDefined {
      public UnaryQOp() {
        Defined<bool, D>.Add(this);
      }
    }

    public static readonly UnaryOp Abs = new UnaryQOp() {
      Id = "Abs",
      MonoOp = NumericOps<T, D1>.Abs,
      Op = NumericOps<double, D1>.Abs,
      CustomExpand = (arg0) => arg0.LessThan(default(T).Bl<T, D>()).Condition(arg0.Negate(), arg0),
    };
    public static readonly UnaryOp Negate = new UnaryPOp() {
      Id = "-",
      MonoOp = NumericOps<T, D1>.Negate,
      Op = NumericOps<double, D1>.Negate,
      CustomInverse0 = (arg0, rhs) => rhs.Negate(),
    };
    public static readonly UnaryOp Recip = typeof(T) == typeof(double) ? new UnaryQOp() {
      Id = "1/",
      MonoOp = NumericOps<T, D1>.Recip,
      Op = NumericOps<double, D1>.Recip,
      CustomInverse0 = (arg0, rhs) => rhs.Recip(),
    } : null;

    public static readonly BinaryOp Add = new BinaryOp() {
      Id = "+",
      MonoOp = NumericOps<T, D1>.Add,
      Op = NumericOps<double, D1>.Add,
      CustomInverse0 = (arg0, arg1, rhs) => rhs.Subtract(arg1),
      CustomInverse1 = (arg0, arg1, rhs) => rhs.Subtract(arg0),
    };
    public static readonly BinaryOp Subtract = new BinaryOp() {
      Id = "-",
      MonoOp = NumericOps<T, D1>.Subtract,
      Op = NumericOps<double, D1>.Subtract,
      CustomInverse0 = (arg0, arg1, rhs) => rhs.Add(arg1),
      CustomInverse1 = (arg0, arg1, rhs) => rhs.Add(arg0).Negate(),
    };

    public static readonly BinaryOp Multiply = new BinaryOp() {
      Id = "*",
      MonoOp = NumericOps<T, D1>.Multiply,
      Op = NumericOps<double, D1>.Multiply,
      CustomInverse0 = (arg0, arg1, rhs) => rhs.Divide(arg1),
      CustomInverse1 = (arg0, arg1, rhs) => rhs.Divide(arg0),
    };
    public static readonly BinaryOp Divide = new BinaryOp() {
      Id = "/",
      MonoOp = NumericOps<T, D1>.Divide,
      Op = NumericOps<double, D1>.Divide,
      CustomInverse0 = (arg0, arg1, rhs) => rhs.Multiply(arg1),
      CustomInverse1 = (arg0, arg1, rhs) => rhs.Divide(arg0).Recip(),
    };
    public static readonly BinaryOp Mod = new BinaryOp() {
      Id = "%",
      MonoOp = NumericOps<T, D1>.Mod,
      Op = NumericOps<double, D1>.Mod,
    };

    public static readonly BinaryOp Min = new BinaryOp() {
      Id = "Min",
      MonoOp = NumericOps<T, D1>.Min,
      Op = NumericOps<double, D1>.Min,
      CustomExpandIfOne = (arg0, arg1) => arg0.LessThan(arg1).Condition(arg0, arg1),
    };
    public static readonly BinaryOp Max = new BinaryOp() {
      Id = "Max",
      MonoOp = NumericOps<T, D1>.Max,
      Op = NumericOps<double, D1>.Max,
      CustomExpandIfOne = (arg0, arg1) => arg0.LessThan(arg1).Condition(arg1, arg0),
    };
    public static readonly TrinaryOp Clamp = new TrinaryOp() {
      Id = "Clamp",
      MonoOp = NumericOps<T, D1>.Clamp,
      Op = NumericOps<double, D1>.Clamp,
      CustomExpandIfOne = (arg0, arg1, arg2) =>
        arg0.Max(arg1).Min(arg2),
    };

  }
  public static partial class DoubleOps<D> where D : Dim<D> {
    public class BinaryOp : Binary2Op<double, double, double, D> {
      public override IOperator Op {
        get {
          return (base.Op == null ? MonoOp : base.Op);
        }
        set {
          base.Op = value;
        }
      }
    }
    public class UnaryOp : UnaryOp<double, double, D>, IShaderDefined, IMathDefined {
      public UnaryOp() {
        Defined<double,D>.Add(this);
      }
      public override IOperator Op {
        get {
          return (base.Op == null ? MonoOp : base.Op);
        }
        set {
          base.Op = value;
        }
      }
    }

    public class Binary2Op : BinaryOp<IVec<double, D>, IVec<double, D>, IVec<double, D1>>, IShaderDefined {
      public Binary2Op() {
        Defined<double,D>.Add(this);
      }
    }
    public class Unary2Op : UnaryOp<IVec<double, D>, IVec<double, D1>>, IShaderDefined, IMathDefined {
      public Unary2Op() {
        Defined<double, D>.Add(this);
      }
    }
    public class Unary3Op : UnaryOp<IVec<double, D>, IVec<bool, D1>>, IShaderDefined, IMathDefined {
      public Unary3Op() {
        Defined<double, D>.Add(this);
      }
      public override IOperator Op {
        get {
          return (base.Op == null ? MonoOp : base.Op);
        }
        set {
          base.Op = value;
        }
      }
      public Func<Expr<IVec<double, D>>, Expr<IVec<bool, D1>>> CustomExpandIfOne {
        set {
          if (Dim<D>.Count == 1) CustomExpand = value;
        }
      }

      public DoubleOps<D1>.Unary3Op MonoOp { private get; set; }
      public override Expr<IVec<bool, D1>> Expand(Expr<IVec<double, D>> Arg0) {
        var e = base.Expand(Arg0);
        if (e != null || Dim<D>.Count == 1 || MonoOp == null) return e;
        var result = false.Repeat<bool>();
        for (int i = 0; i < Dim<D>.Count; i++) 
          result = result.Or(MonoOp.Make(Arg0.Access(0).Bl0()));
        return result;
      }
    }
    public class TrinaryOp : TrinaryOp<double, double, double, double, D>, IShaderDefined {
      public TrinaryOp() {
        Defined<double,D>.Add(this);
      }
      public override IOperator Op {
        get {
          return (base.Op == null ? MonoOp : base.Op);
        }
        set {
          base.Op = value;
        }
      }
    }
    public class Trinary2Op : TrinaryOp<IVec<double,D>,IVec<double,D>,IVec<double,D1>,IVec<double,D>>, IShaderDefined {
      public Trinary2Op() {
        Defined<double,D>.Add(this);
      }
    }

    public static readonly UnaryOp Saturate = new UnaryOp() {
      Id = "Saturate",
      MonoOp = DoubleOps<D1>.Saturate,
      CustomExpand = (arg0) =>
        arg0.Clamp(0d.Bl<double, D>(), 1d.Bl<double, D>()),
    };

    public static readonly UnaryOp Floor = new UnaryOp() {
      Id = "Floor",
      MonoOp = DoubleOps<D1>.Floor,
      CustomExpand = (arg0) => arg0.Subtract(arg0.Mod(1d.Bl<double, D>())),
    };
    public static readonly UnaryOp Truncate = new UnaryOp() {
      Id = "Truncate", Alt= "trunc",
      MonoOp = DoubleOps<D1>.Truncate,
      // XXX: how are floor and truncate different?
      CustomExpand = (arg0) => arg0.Subtract(arg0.Mod(1d.Bl<double, D>())),
    };
    public static readonly UnaryOp Ceiling = new UnaryOp() {
      Id = "Ceiling", Alt = "ceil",
      MonoOp = DoubleOps<D1>.Ceiling,
      CustomExpand = (arg0) => arg0.Add(1d.Bl<double, D>().Subtract(arg0.Mod(1d.Bl<double, D>()))),
    };
    public static readonly UnaryOp Round = new UnaryOp() {
      Id = "Round",
      MonoOp = DoubleOps<D1>.Round,
      CustomExpand = (arg0) =>
        arg0.Multiply(2d.Bl<double, D>()).Floor().Divide(2d.Bl<double, D>()).Ceiling(),
    };
    public static readonly UnaryOp Frac = new UnaryOp() {
      Id = "Frac",
      MonoOp = DoubleOps<D1>.Frac,
      CustomExpand = (arg0) =>
        arg0.Mod(1d.Bl<double, D>()).Abs(),
    };
    public static readonly BinaryOp Fmod = new BinaryOp() {
      Id = "Fmod",
      MonoOp = DoubleOps<D1>.Fmod,
    };
    public static readonly BinaryOp FrExp = new BinaryOp() {
      Id = "FrExp",
      MonoOp = DoubleOps<D1>.FrExp,
    };
    public static readonly UnaryOp Square = new UnaryOp() {
      Id = "Square",
      MonoOp = DoubleOps<D1>.Square,
      CustomExpand = (arg0) => arg0.Multiply(arg0),
      CustomInverse0 = (arg0, rhs) => rhs.Sqrt(),
    };
    public static readonly UnaryOp Sqrt = new UnaryOp() {
      Id = "Sqrt",
      MonoOp = DoubleOps<D1>.Sqrt,
      CustomInverse0 = (arg0, rhs) => rhs.Square(),
    };
    public static readonly UnaryOp Rsqrt = new UnaryOp() {
      Id = "Rsqrt",
      MonoOp = DoubleOps<D1>.Rsqrt,
      CustomExpand = (arg0) => arg0.Sqrt().Recip(),
      CustomInverse0 = (arg0, rhs) => rhs.Square().Recip(),
    };

    public static readonly Binary2Op Dot = new Binary2Op() {
      Id = "Dot",
      Op = DoubleOps<D1>.Dot,
      CustomExpand = (arg0, arg1) => {
        var result = 0d.Repeat<double>();
        for (int i = 0; i < Dim<D>.Count; i++)
          result = result.Add(arg0.Access(i).Bl0().Multiply(arg1.Access(i).Bl0()));
        return result;
      },
    };
    public static readonly Unary2Op LengthSquared = new Unary2Op() {
      Id = "LengthSquared",
      Op = DoubleOps<D1>.LengthSquared,
      CustomExpand = (arg0) => arg0.Dot(arg0),
    };
    public static readonly Unary2Op Length = new Unary2Op() {
      Id = "Length",
      Op = DoubleOps<D1>.Length,
      CustomExpand = (arg0) => arg0.LengthSquared().Sqrt(),
    };
    public static readonly UnaryOp Normalize = new UnaryOp() {
      Id = "Normalize",
      MonoOp = DoubleOps<D1>.Normalize,
      CustomExpand = (arg0) => arg0.Divide(arg0.Length().Bl<double, D>()),
    };
    public static readonly Binary2Op Distance = new Binary2Op() {
      Id = "Distance",
      Op = DoubleOps<D1>.Distance,
      CustomExpand = (arg0, arg1) => arg0.Subtract(arg1).Length(),
    };

    public static readonly UnaryOp Ddx = new UnaryOp() {
      Id = "Ddx",
      MonoOp = DoubleOps<D1>.Ddx,
    };
    public static readonly UnaryOp Ddy = new UnaryOp() {
      Id = "Ddy",
      MonoOp = DoubleOps<D1>.Ddy,
    };
    public static readonly UnaryOp FWidth = new UnaryOp() {
      Id = "FWidth",
      MonoOp = DoubleOps<D1>.FWidth,
      CustomExpand = (arg0) => arg0.Ddx().Abs().Add(arg0.Ddy().Abs()),
    };

    public static readonly UnaryOp Exp = new UnaryOp() {
      Id = "Exp",
      MonoOp = DoubleOps<D1>.Exp,
      CustomExpand = (arg0) => Math.E.Bl<double, D>().Pow(arg0),
      CustomInverse0 = (arg0, rhs) => rhs.Log(),
    };
    public static readonly UnaryOp Exp2 = new UnaryOp() {
      Id = "Exp2",
      MonoOp = DoubleOps<D1>.Exp2,
      CustomExpand = (arg0) => 2d.Bl<double, D>().Pow(arg0),
      CustomInverse0 = (arg0, rhs) => rhs.Log2(),
    };
    public static readonly UnaryOp Exp10 = new UnaryOp() {
      Id = "Exp10",
      MonoOp = DoubleOps<D1>.Exp10,
      CustomExpand = (arg0) => 10d.Bl<double, D>().Pow(arg0),
      CustomInverse0 = (arg0, rhs) => rhs.Log10(),
    };
    public static readonly BinaryOp Pow = new BinaryOp() {
      Id = "Pow",
      MonoOp = DoubleOps<D1>.Pow,
      CustomInverse1 = (arg0, arg1, rhs) => rhs.LogN(arg0),
    };

    public static readonly UnaryOp Log = new UnaryOp() {
      Id = "Log",
      MonoOp = DoubleOps<D1>.Log,
      CustomExpand = (arg0) => arg0.LogN(Math.E.Bl<double, D>()),
      CustomInverse0 = (arg0, rhs) => rhs.Exp(),
    };
    public static readonly UnaryOp Log2 = new UnaryOp() {
      Id = "Log2",
      MonoOp = DoubleOps<D1>.Log2,
      CustomExpand = (arg0) => arg0.LogN(2d.Bl<double, D>()),
      CustomInverse0 = (arg0, rhs) => rhs.Exp2(),
    };
    public static readonly UnaryOp Log10 = new UnaryOp() {
      Id = "Log10",
      MonoOp = DoubleOps<D1>.Log10,
      CustomExpand = (arg0) => arg0.LogN(10d.Bl<double, D>()),
      CustomInverse0 = (arg0, rhs) => rhs.Exp10(),
    };
    public static readonly BinaryOp LogN = new BinaryOp() {
      Id = "Log",
      MonoOp = DoubleOps<D1>.LogN,
      CustomInverse0 = (arg0, arg1, rhs) => arg1.Pow(rhs),
    };
    public static readonly BinaryOp Ldexp = new BinaryOp() {
      Id = "Ldexp",
      MonoOp = DoubleOps<D1>.Ldexp,
      CustomExpand = (arg0, arg1) => arg0.Multiply(arg1.Exp2()),
    };
    public static readonly UnaryOp Sign = new UnaryOp() {
      Id = "Sign",
      MonoOp = DoubleOps<D1>.Sign,
      CustomExpand = (arg0) => arg0.LessThan(0d.Bl<double,D>()).Condition<double,D>((-1d).Bl<double,D>(), (+1d).Bl<double,D>()),
    };

    public static readonly TrinaryOp Faceforward = new TrinaryOp() {
      Id = "Faceforward",
      MonoOp = DoubleOps<D1>.Faceforward,
      //−n × sign(dot(i, ng)).
      CustomExpand = (n, i, ng) => n.Negate().Multiply(i.Dot(ng).Sign().Bl<double,D>()),
    };
    public static readonly BinaryOp Reflect = new BinaryOp() {
      Id = "Reflect",
      MonoOp = DoubleOps<D1>.Reflect,
      // v = i - 2 * n * dot(i•n) .
      CustomExpand = (i, n) => i.Subtract(2d.Bl<double,D>().Multiply(n).Multiply(i.Dot(n).Bl<double,D>())),
    };
    public static readonly Trinary2Op Refract = new Trinary2Op() {
      Id = "Refract",
      /*
k = 1.0 - eta * eta * (1.0 - dot(N, I) * dot(N, I))
if (k < 0.0)
return genType(0.0)
else
return eta * I - (eta * dot(N, I) + sqrt(k)) * N
       */
      CustomExpand = (i, n, eta) => {
        var one = 1.0.Repeat<double>();
        var zero = 0.0.Repeat<double>();
        var k = eta.Square().Multiply(one.Subtract(n.Dot(i)));
        k = one.Subtract(k);

        var k2 = (eta.Multiply(n.Dot(i).Add(k.Sqrt())).Bl<double, D>().Multiply(n));
        k2 = eta.Bl<double, D>().Multiply(i).Subtract(k2);

        return (k.LessThan(zero)).Bl<bool,D>().Condition(zero.Bl<double, D>(), k2);
      },
      Op = DoubleOps<D1>.Refract,
    };

    public static readonly Unary3Op IsNaN = new Unary3Op() {
      MonoOp = DoubleOps<D1>.IsNaN,
      Id = "IsNaN",
      CustomExpandIfOne = (arg0) => arg0.Equals0(double.NaN.Bl<double,D>()).Any(),
    };
    public static readonly Unary3Op IsInfinity = new Unary3Op() {
      MonoOp = DoubleOps<D1>.IsInfinity,
      Id = "IsInfinity", Alt = "isinf",
      CustomExpandIfOne = (arg0) => arg0.IsPositiveInfinity().Or(arg0.IsNegativeInfinity()),
    };
    public static readonly Unary3Op IsFinite = new Unary3Op() {
      MonoOp = DoubleOps<D1>.IsFinite,
      Id = "IsFinite", Alt = "isfin",
      CustomExpandIfOne = (arg0) => arg0.IsInfinity().Not(),
    };
    public static readonly Unary3Op IsPositiveInfinity = new Unary3Op() {
      MonoOp = DoubleOps<D1>.IsPositiveInfinity,
      Id = "IsPositiveInfinity",
      CustomExpandIfOne = (arg0) => arg0.Equals0(double.PositiveInfinity.Bl<double, D>()).Any(),
    };
    public static readonly Unary3Op IsNegativeInfinity = new Unary3Op() {
      MonoOp = DoubleOps<D1>.IsNegativeInfinity,
      Id = "IsNegativeInfinity",
      CustomExpandIfOne = (arg0) => arg0.Equals0(double.NegativeInfinity.Bl<double, D>()).Any(),
    };

  }
  public static partial class TrigOps<D> where D : Dim<D> {
    public class BinaryOp : DoubleOps<D>.BinaryOp {}
    public class UnaryOp : DoubleOps<D>.UnaryOp { }
    public class Binary2Op : DoubleOps<D>.Binary2Op { }
    public class Unary2Op : DoubleOps<D>.Unary2Op { }
    public class TrinaryOp : DoubleOps<D>.TrinaryOp { }

    public static readonly UnaryOp Sin = new UnaryOp() {
      MonoOp = TrigOps<D1>.Sin, 
      Id = "Sin",
      CustomInverse0 = (arg, result) => result.Asin(),
    };
    public static readonly UnaryOp Asin = new UnaryOp() {
      MonoOp = TrigOps<D1>.Asin,
      Id = "Asin",
      CustomInverse0 = (arg, result) => result.Sin(),
    };
    public static readonly UnaryOp Sinh = new UnaryOp() {
      MonoOp = TrigOps<D1>.Sinh,
      Id = "Sinh",
    };
    public static readonly UnaryOp Cos = new UnaryOp() {
      MonoOp = TrigOps<D1>.Cos,
      Id = "Cos",
      CustomInverse0 = (arg, result) => result.Acos(),
    };
    public static readonly UnaryOp Cosh = new UnaryOp() {
      MonoOp = TrigOps<D1>.Cosh,
      Id = "Cosh",
    };
    public static readonly UnaryOp Acos = new UnaryOp() {
      MonoOp = TrigOps<D1>.Acos,
      Id = "Acos",
      CustomInverse0 = (arg, result) => result.Cos(),
    };
    public static readonly UnaryOp Tan = new UnaryOp() {
      MonoOp = TrigOps<D1>.Tan,
      Id = "Tan",
      CustomInverse0 = (arg, result) => result.Atan(),
    };
    public static readonly UnaryOp Atan = new UnaryOp() {
      MonoOp = TrigOps<D1>.Atan,
      Id = "Atan",
      CustomInverse0 = (arg,result) => result.Tan(),
    };
    public static readonly BinaryOp Atan2 = new BinaryOp() {
      MonoOp = TrigOps<D1>.Atan2,
      Id = "Atan2",
    };
    public static readonly UnaryOp Radians = new UnaryOp() {
      MonoOp = TrigOps<D1>.Radians,
      Id = "Radians",
      CustomExpand = (Degrees) =>
        Degrees.Divide((180.0 / Math.PI).Bl<double,D>()),
      CustomInverse0 = (arg,Radians) => Radians.Degrees(),
    };
    public static readonly UnaryOp Degrees = new UnaryOp() {
      MonoOp = TrigOps<D1>.Degrees,
      Id = "Degrees",
      CustomExpand = (Radians) =>
        Radians.Multiply((180.0 / Math.PI).Bl<double, D>()),
      CustomInverse0 = (arg,Degrees) => Degrees.Radians(),
    };
  }
  public static partial class MatrixOps<D>
    where D  : Dim2D<D> {
    public static int IndexFor(int m, int n) {
      return m * Dim2D<D>.NCount + n;
    }
    public static Expr<int> IndexFor(Expr<int> m, Expr<int> n) {
      return (m.Bl() * Dim2D<D>.NCount + n.Bl()).Underlying;
    }
  }
  public static partial class MatrixOps<T, D>
    where D : Dim2D<D> {
    public class UnaryOp : UnaryOp<T,T,D> {
    }
    public static readonly UnaryOp Transpose = new UnaryOp() {
      Id = "Transpose",
      CustomInverse0 = (arg, result) => result.Transpose<T, D>(),
      CustomExpand = (arg0) => {
        var Es = new Expr<T>[Dim2D<D>.Count];
        for (int i = 0; i < Dim2D<D>.MCount; i++) {
          for (int j = 0; j < Dim2D<D>.NCount; j++) {
            Es[MatrixOps<D>.IndexFor(j, i)] = arg0.Access<T, D>(i, j);
          }
        }
        return Composite<T, D>.Make(Es);
      },
      Op = MatrixOps<double,Dim2D<D4,D4>>.Transpose,
    };
  }
  public static partial class MatrixOps<T, M, N, P>
    where M : Dim1D<M>
    where N : Dim1D<N>
    where P : Dim1D<P>  {
    static MatrixOps() {
    }
    public class MultOp : BinaryOp<IVec<T, Dim2D<M, N>>, IVec<T, Dim2D<N, P>>, IVec<T, Dim2D<M, P>>> { }
    public static readonly MultOp Mul = new MultOp() {
      Id = "Mul",
      CustomExpand = (mn, np) => {
        var Es = new Expr<T>[Dim<Dim2D<M,P>>.Count];
        for (int m = 0; m < Dim<M>.Count; m++) {
          for (int p = 0; p < Dim<P>.Count; p++) {
            var result = default(T).Repeat<T>();
            for (int n = 0; n < Dim<N>.Count; n++) {
              var mn0 = mn.Access(m, n).Bl0<T>();
              var np0 = np.Access(n, p).Bl0<T>();
              result = result.Add(mn0.Multiply(np0));
            }
            Es[MatrixOps<Dim2D<M,P>>.IndexFor(m, p)] = result.Access(0);
          }
        }
        return Composite<T, Dim2D<M, P>>.Make(Es);
      },
      Op = MatrixOps<double,D4,D4,D4>.Mul,
    };
  }

  public static partial class SquareMatrixOps<S>
    where S : Dim1D<S> {

    public class Unary2Op : DoubleOps<Dim2D<S,S>>.Unary2Op { }

    public class TransformOp : BinaryOp<IVec<double, Dim2D<S, S>>, IVec<double, S>, IVec<double, S>> {}


    /*
    public static readonly TransformOp Transform = new TransformOp() {
      CustomExpand = (matrix, point) => {
        var Es = new Expr<double>[Dim<S>.Count];
        for (int i = 0; i < Es.Length; i++) {
        }
      },
    };
    public static BRAND Transform<T, BRAND, ARITY>(this MatrixBl<D4, D4> M, Point3DBl<T, BRAND, ARITY> P)
      where BRAND : Point3DBl<T, BRAND, ARITY>
      where ARITY : Arity3<double, T, ARITY>, new() {

      var newX = P.X * M.M00() + P.Y * M.M10() + P.Z * M.M20() + M.OffsetX();
      var newY = P.X * M.M01() + P.Y * M.M11() + P.Z * M.M21() + M.OffsetY();
      var newZ = P.X * M.M02() + P.Y * M.M12() + P.Z * M.M22() + M.OffsetZ();
      var ret = Point3DBl<T, BRAND, ARITY>.Make(newX, newY, newZ);
      var w = P.X * M.M03() + P.Y * M.M13() + P.Z * M.M23() + M.M33();
      return ret / w;
    }
*/
    public static readonly Unary2Op Determinant = new Unary2Op() {
      Id = "Determinant",
      CustomExpand = (matrix) => {
        if (Dim<S>.Count == 1) return matrix.Access(0, 0).Bl0<double>();
        else if (Dim<S>.Count == 2) {
          var a = matrix.Access(0, 0).Bl0<double>();
          var b = matrix.Access(0, 1).Bl0<double>();
          var c = matrix.Access(1, 0).Bl0<double>();
          var d = matrix.Access(1, 1).Bl0<double>();
          return a.Multiply(d).Subtract(b.Multiply(c));
        } else if (Dim<S>.Count == 3) {
          var a = matrix.Access(0, 0).Bl0<double>();
          var b = matrix.Access(0, 1).Bl0<double>();
          var c = matrix.Access(0, 2).Bl0<double>();
          var d = matrix.Access(1, 0).Bl0<double>();
          var e = matrix.Access(1, 1).Bl0<double>();
          var f = matrix.Access(1, 2).Bl0<double>();
          var g = matrix.Access(2, 0).Bl0<double>();
          var h = matrix.Access(2, 1).Bl0<double>();
          var i = matrix.Access(2, 2).Bl0<double>();
          //det (A) = aei − afh − bdi + bfg + cdh − ceg

          var aei = a.Multiply(e).Multiply(i);
          var afh = a.Multiply(f).Multiply(h);
          var bdi = b.Multiply(d).Multiply(i);
          var bfg = b.Multiply(f).Multiply(g);
          var cdh = c.Multiply(d).Multiply(h);
          var ceg = c.Multiply(e).Multiply(g);
          return aei.Subtract(afh).Subtract(bdi).Add(bfg).Add(cdh).Subtract(ceg);
        } else if (Dim<S>.Count == 4) {
//det(A) = a11a22a33a44 - a11a22a34a43 + a11a23a34a42 - a11a23a32a44 + a11a24a32a43 - a11a24a33a42
  //     - a12a23a34a41 + a12a23a31a44 - a12a24a31a43 + a12a24a33a41 - a12a21a33a44 + a12a21a34a43
    //   + a13a24a31a42 - a13a24a32a41 + a13a21a32a44 - a13a21a34a42 + a13a22a34a41 - a13a22a31a44
      // - a14a21a32a43 + a14a21a33a42 - a14a22a33a41 + a14a22a31a43 - a14a23a31a42 + a14a23a32a41

          var result = 0d.Repeat<double>();
          for (int i = 0; i < 4; i++) {
            var signI = (i % 2 == 0 ? +1d : -1d);
            for (int j = 0; j < 6; j++) {
              var signJ = (j % 2 == 0 ? +1d : -1d);
              var sign = (signI * signJ).Repeat<double>();
              var t0 = matrix.Access(0, (0 + i + (new int[] { 0, 0, 0, 0, 0, 0 })[j]) % 4).Bl0<double>();
              var t1 = matrix.Access(1, (1 + i + (new int[] { 0, 0, 1, 1, 2, 2 })[j]) % 4).Bl0<double>();
              var t2 = matrix.Access(2, (2 + i + (new int[] { 0, 1, 1, 3, 3, 4 })[j]) % 4).Bl0<double>();
              var t3 = matrix.Access(3, (3 + i + (new int[] { 0, 3, 2, 0, 3, 2 })[j]) % 4).Bl0<double>();
              result = result.Add(sign.Multiply(t0).Multiply(t1).Multiply(t2).Multiply(t3));
            }
          }
          return result;
        }
        throw new NotImplementedException();
      },
    };
  }
}
